function gccolor(map,lon,lat,data,varargin)

% _________________________________________________________________________
% Call:
% gccolor(map,lon,lat,data)
%
% Input:
%   Needed:
%   - map : 'cl' for a classical pcolor
%           'clc' for a classical contourf
%           'mp' for a m_map plot (pcolor)
%           'mpc' for a m_map plot (contourf)
%           'mquiver' for a quiver map (m_quiver)
%           'quiver' for a quiver plot (quiver)
%   - lon : x coordinates
%   - lat : y coordinates
%   - data : Data to plot (for quiver: data <=> [U V])
%
%   List of options:
%   - 'figure' : 'new' (default) /'hold'
%       -> new: open a new figure for the map
%       -> hold: use the current figure for the map
%   - 'coastal' : 'on'/'GdG'/'off' (default) 
%       -> on = high resolution coast line (slower computation)
%       -> GdG = high resolution coast line of Bay of Biscay
%       -> off = low resolution coast line
%   - 'clb' : 'on' (default) /'off' -> on = colorbar, off = no colorbar
%   - 'caxis' : [mincaxis maxcaxis]
%       -> mincaxis : Lower limit for the caxis
%       -> maxcaxis : Higer limit for the caxis
%   - 'title' -> Title to add (string)
%   - 'domain' : [milat malat milon malon]
%       -> milat : Min Latitude
%       -> malat : Max Latitude
%       -> milon : Min Longitude
%       -> malon : Max Longitude
%   - 'legendlocation' : - ONLY FOR QUIVER PLOT - to specify the reference 
%                        vector location
%       -> 1: bottom/left corner
%       -> 2: bottom/right corner
%       -> 3: top/right corner
%       -> 4: top/left corner
%   - 'contours' : in the case of contour plot, you can define contours
%   values (it does not affect other kind of plots)
%   - 'visible' : 'on' (default) /'off' -> figure visible on the screen 
%                                          or not
%   - 'centeredpixel' : 'off' (default) /'on' (cl and mp map)
%       -> By default, using pcolor, matlab takes coordinates as the
%       bottom/left corner. When this option is switched to 'on', 
%       coordinates are used as the center of the pixel.
%       (Only for 1D or regular grids)
%
%   Example:
%   % -- Map with a high resolution coastal line
%   gccolor('mp',lon,lat,data,'coastal','on') 
%
% Created: G. Charria, 12/2007
% Last update: 02/2010
% -----------------------------------

if nargin < 1, help(mfilename), return, end
narg=nargin - 4;
if mod(narg,2) ~=0 
    error('Wrong number of input arguments (see help gccolor) !')   
end

switch map
    case {'cl','clc','mp','mpc'}
        % Add a NaN line and a NaN column
        data(:,size(data,2)+1)=NaN;
        data(size(data,1)+1,:)=NaN;
        %size(data)
        if ((size(lon,1)>1) && (size(lon,2)>1))
            lon(size(lon,1)+1,:)=NaN;
            %size(lon)
            lon(:,size(lon,2)+1)=lon(:,size(lon,2))+(lon(:,size(lon,2))-lon(:,size(lon,2)-1));
            lon(size(lon,1),:)=lon(size(lon,1)-1,:)+(lon(size(lon,1)-1,:)-lon(size(lon,1)-2,:));
            %size(lon)
            lat(size(lat,1)+1,:)=NaN;
            lat(:,size(lat,2)+1)=lat(:,size(lat,2))+(lat(:,size(lat,2))-lat(:,size(lat,2)-1));
            lat(size(lat,1),:)=lat(size(lat,1)-1,:)+(lat(size(lat,1)-1,:)-lat(size(lat,1)-2,:));
            %size(lat)
            % - Test the data dimension compared to lon,lat dimensions
            if (size(lon,1)==size(data,2))
                data=data';
            end
        elseif (size(lon,1))
            lon(length(lon)+1)=lon(length(lon))+(lon(length(lon))-lon(length(lon)-1));
            lat(length(lat)+1)=lat(length(lat))+(lat(length(lat))-lat(length(lat)-1));
        elseif (size(lon,2))
            lon(length(lon)+1)=lon(length(lon))+(lon(length(lon))-lon(length(lon)-1));
            lat(length(lat)+1)=lat(length(lat))+(lat(length(lat))-lat(length(lat)-1));
            data=data';
        else
            error('Error in lon and/or lat size !')
        end
        % - Default parameters
        clb='on';
    case {'quiver','mquiver'}
        if size(lon)==size(lat) & size(lon,1)==size(data,1)
            U=data(:,1:size(lon,2));
            V=data(:,size(lon,2)+1:end);
        else
            error('Error in the input variable dimensions !')
        end 
        % - Default parameters
        clb='off'; 
end

% -- Scan options

% - Default parameters
caxison=[];
titre=[];
domain=[];
coastal='off';
figtest='new';
legendlocation=1;
ctrs=[];
visible='on';
centeredpixel='off';

for iopt = 1 : 2 : narg
    optname  = varargin{iopt};
    optvalue = varargin{iopt+1};
    switch lower(optname)
        case 'clb'
            clb=optvalue;
        case 'caxis'
            caxison=optvalue;
        case 'title'
            titre=optvalue;
        case 'domain'
            domain=optvalue;
        case 'coastal'
            coastal=optvalue;
        case 'figure'
            figtest=optvalue;
        case 'title'
            titre=optvalue;
        case 'legendlocation'
            legendlocation=optvalue;
        case 'contours'
            ctrs=optvalue;
        case 'visible'
            visible=optvalue;
        case 'centeredpixel'
            centeredpixel='on';
        otherwise
            error('Wrong option, please check help gccolor.')
            
    end
end


switch figtest
    case 'new'
        figure   
end
set(gcf,'Visible',visible)
switch map
    case 'mp'
        if strcmp('on',centeredpixel) % centered pixels
            if ((size(lon,1)>1) && (size(lon,2)>1)) 
                
                
                if max(max(diff(lon,1,1)))==min(min(diff(lon,1,1))) &&...
                        max(max(diff(lon,1,2)))==min(min(diff(lon,1,2))) && ...
                        max(max(diff(lat,1,1)))==min(min(diff(lat,1,1))) && ...
                        max(max(diff(lat,1,2)))==min(min(diff(lat,1,2)))                                        
                    
                    xp1=lon(2:end,2:end);
                    x0=lon(1:end-1,1:end-1);
                    lon(1:end-1,1:end-1)=x0-((xp1-x0)/2);
                    yp1=lat(2:end,2:end);
                    y0=lat(1:end-1,1:end-1);
                    lat(1:end-1,1:end-1)=y0-((yp1-y0)/2);
                    
                    
                    lat(end,:)=lat(end-1,:)+(lat(end-1,:)-lat(end-2,:));
                    lat(:,end)=lat(:,end-1)+(lat(:,end-1)-lat(:,end-2));
                    lon(end,:)=lon(end-1,:)+(lon(end-1,:)-lon(end-2,:));
                    lon(:,end)=lon(:,end-1)+(lon(:,end-1)-lon(:,end-2));
                    
                else
                    
                    % Under dev in test_grid_pcolor.m
                    disp('This option do not work with 2D irregular grid !')
                    return
                end
                
                
            else
                dlon=diff(lon);
                lon(2:end)=lon(2:end)-(dlon/2);
                lon(1)=lon(1)-(dlon(1)/2);
                dlat=diff(lat);
                lat(2:end)=lat(2:end)-(dlat/2);
                lat(1)=lat(1)-(dlat(1)/2);  

            end
        end
        if ~isempty(domain)
            m_proj('equidistant','lat',[domain(1) domain(2)],'lon',[domain(3) domain(4)])
        else
            m_proj('equidistant','lat',[min(min(lat)) max(max(lat))],'lon',[min(min(lon)) max(max(lon))]);
        end
        m_pcolor(lon,lat,data);
        shading flat
        switch coastal
            case 'GdG'
                if ~exist('TCote.mat','file'), m_gshhs_f('save','TCote');end;
                m_usercoast('TCote','patch',[.6 .6 .6],'edgecolor','k');
            case 'on'
                m_gshhs_f('patch',[.6 .6 .6]);
            otherwise
                m_coast('patch',[.6 .6 .6]);
        end            
        m_grid('linestyle',':','box','fancy','tickdir','in','Fontsize',12);
    case 'mpc'
        if ~isempty(domain)
            m_proj('equidistant','lat',[domain(1) domain(2)],'lon',[domain(3) domain(4)])
        else
            m_proj('equidistant','lat',[min(min(lat)) max(max(lat))],'lon',[min(min(lon)) max(max(lon))]);
        end
        m_contourf(lon,lat,data,64);
        shading flat
        switch coastal
            case 'GdG'
                if ~exist('TCote.mat','file'), m_gshhs_f('save','TCote');end;
                m_usercoast('TCote','patch',[.6 .6 .6],'edgecolor','k');
            case 'on'
                m_gshhs_f('patch',[.6 .6 .6]);
            otherwise
                m_coast('patch',[.6 .6 .6]);
        end  
        m_grid('linestyle',':','box','fancy','tickdir','in','Fontsize',12);
    case 'cl'
        if strcmp('on',centeredpixel)
            if ((size(lon,1)>1) && (size(lon,2)>1))
                
                if max(max(diff(lon,1,1)))==min(min(diff(lon,1,1))) &&...
                        max(max(diff(lon,1,2)))==min(min(diff(lon,1,2))) && ...
                        max(max(diff(lat,1,1)))==min(min(diff(lat,1,1))) && ...
                        max(max(diff(lat,1,2)))==min(min(diff(lat,1,2)))                                        
                    
                    xp1=lon(2:end,2:end);
                    x0=lon(1:end-1,1:end-1);
                    lon(1:end-1,1:end-1)=x0-((xp1-x0)/2);
                    yp1=lat(2:end,2:end);
                    y0=lat(1:end-1,1:end-1);
                    lat(1:end-1,1:end-1)=y0-((yp1-y0)/2);
                    
                    
                    lat(end,:)=lat(end-1,:)+(lat(end-1,:)-lat(end-2,:));
                    lat(:,end)=lat(:,end-1)+(lat(:,end-1)-lat(:,end-2));
                    lon(end,:)=lon(end-1,:)+(lon(end-1,:)-lon(end-2,:));
                    lon(:,end)=lon(:,end-1)+(lon(:,end-1)-lon(:,end-2));
                    
                else
                    
                    % Under dev in test_grid_pcolor.m
                    disp('This option do not work with 2D irregular grid !')
                    return
                end
                
            else
                dlon=diff(lon);
                lon(2:end)=lon(2:end)-(dlon/2);
                lon(1)=lon(1)-(dlon(1)/2);
                dlat=diff(lat);
                lat(2:end)=lat(2:end)-(dlat/2);
                lat(1)=lat(1)-(dlat(1)/2);  

            end
        end
        pcolor(lon,lat,data)
        if ~isempty(domain)
            set(gca,'YLim',[domain(1) domain(2)],'XLim',[domain(3) domain(4)],'Fontsize',12)
        end
        shading flat
    case 'clc'
        if ~isempty(ctrs)
            contourf(lon,lat,data,ctrs)
        else
            contourf(lon,lat,data,64)
        end
        if ~isempty(domain)
            set(gca,'YLim',[domain(1) domain(2)],'XLim',[domain(3) domain(4)],'Fontsize',12)
        end
        shading flat
    case 'mquiver'
        if ~isempty(domain)
            m_proj('equidistant','lat',[domain(1) domain(2)],'lon',[domain(3) domain(4)])
        else
            m_proj('equidistant','lat',[min(min(lat)) max(max(lat))],'lon',[min(min(lon)) max(max(lon))]);
        end
        switch coastal
            case 'GdG'
                if ~exist('TCote.mat','file'), m_gshhs_f('save','TCote');end;
                m_usercoast('TCote','patch',[.6 .6 .6],'edgecolor','k');
            case 'on'
                m_gshhs_f('patch',[.6 .6 .6]);
            otherwise
                m_coast('patch',[.6 .6 .6]);
        end  
        m_grid('linestyle',':','box','fancy','tickdir','in','Fontsize',12);
        hold on;
        m_quiver(lon,lat,U,V,1); 
%         
%         % legendlocation:
%         % 1: bottom/left corner
%         % 2: bottom/right corner
%         % 3: top/right corner
%         % 4: top/left corner
%         if legendlocation==1
%             % Get the current axis limits
%             xlim=get(gca,'xlim'); xp1=xlim(1); xp2=xlim(2);
%             ylim=get(gca,'ylim'); yp1=ylim(1); yp2=ylim(2);
%             % set padding around the reference vector
%             padx=diff(xlim)/100; 
%             pady=diff(ylim)/100;
%             % Set x position of reference vector
%             xstart=xp1+padx;
%             xend=xstart+scalelength;
%             % Plot reference text in lower right hand corner
%             ht=text(xstart,yp1+pady,reftext,'Visible','off','Parent',gca,'FontSize',8.5,...
%                 'VerticalAlignment','Bottom','HorizontalAlignment','Right');
%             textextent=get(ht,'Extent');            
%             % Draw patch over area of vector key 
%             xl=textextent(1)-padx;
%             xr=xp2;
%             yb=yp1;
%             yt=textextent(2)+textextent(4)+pady;
%             hp=patch([xl; xl; xr; xr],[yb; yt; yt; yb],[2; 2; 2; 2],'w',...
%                 'LineWidth',0.5,'Parent',gca);
%             uistack(hp,'top');
%             % Redraw reference text on top of patch
%             ht=text(xstart,yp1+pady,2.1,reftext,'Parent',gca,'FontSize',8.5,...
%                 'VerticalAlignment','Bottom','HorizontalAlignment','Right');
%             % Set y position of reference vector
%             yend=textextent(2)+textextent(4)/2;
%             ystart=yend;
%         end
%         
%         
        % -- Attention trace vecteur unite --
        %m_quiver(max(max(lon))-(max(max(lon))-min(min(lon)))/4,max(max(lat))-(max(max(lat))-min(min(lon)))/4,10,0,1); 
        %m_vec(1,max(max(lon))-(max(max(lon))-min(min(lon)))/4,max(max(lat))-(max(max(lat))-min(min(lon)))/4,10,0,'key','10 x Unit vector')
    case 'quiver'
        quiver(lon,lat,U,V)
        if ~isempty(domain)
            set(gca,'YLim',[domain(1) domain(2)],'XLim',[domain(3) domain(4)],'Fontsize',12)
        end        
    otherwise
        error('Error in the first function parameter')      
end

if ~isempty(caxison)
    caxis(caxison)
end
if strcmp('on',clb)
    colorbar
end
if ~isempty(titre)
    title(titre,'Fontsize',12)
end
